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Final  Progress  Report 
Problem  studied: 

The  overall  objective  of  this  work  is  to  understand  information  flow  on  a  network  through  the  analogous 
process  of  biological  disease  spread  in  a  noisy  environment.  Currently,  the  PI  is  developing  new  multi-patch 
models  that  address  the  issue  of  scalability,  where  the  size  of  the  groupings  is  determined  by  the  accuracy  and 
scope  of  the  results  needed.  She  has  developed  a  computational  tool  to  predict  changes  in  dynamics  due  to  noise 
probabilistically.  This  method  numerically  approximates  the  mechanisms  of  transport  in  a  mathematical  space, 
and  provides  a  way  to  visualize  it  using  a  matrix  representation.  Commonly,  the  only  data  available  in  the  field 
is  the  number  of  infected  individuals  reported.  Using  embedding  techniques,  this  data  can  emulate  the  dynamics 
of  the  full  system  and  be  used  in  our  methods.  In  practice,  the  PI  has  successfully  compared  a  stochastic 
bifurcation  in  a  controlled  experiment  to  a  theoretical  laser  model.  The  transport  tools  developed  to  date  allow 
control  of  the  system  by  targeting  or  avoiding  high  probability  regions.  The  PI  has  a  working  algorithm 
submitted  for  publication  and  is  planning  its  implementation  in  a  future  experiment. 


Summary  of  the  most  important  results: 

1 .  Development  of  the  continuous  transition  matrix 

One  important  mathematical  question  is  the  precise  effect  of  noise  on  an  ensemble.  We  focus  on  the 
dynamics  of  ensembles  that  obey  deterministic  dynamics,  but  are  subject  to  continuously  added  noise.  Since 
stable  configurations  of  interacting  particles  may  be  thought  of  as  living  in  potential  wells,  changing  the 
dynamics  to  a  new  state  may  require  escape  from  a  given  potential,  and  we  consider  the  probability  of 
escape  happening  due  to  natural  or  imposed  fluctuations.  Escape  is  an  example  of  a  large  fluctuation.  If 
fluctuations  are  small  on  average,  the  ensemble  will  continue  to  emulate  the  current  state.  Simple  time 
simulations  from  models  have  difficulty  generating  rare  events  like  escape  for  proper  analysis.  Therefore, 
we  have  designed  a  new,  efficient  technique  that  can  be  used  approximate  the  probability  density  function 
based  on  spatial  averages.  Our  mathematical  method  is  based  on  a  stochastic  operator  formalism,  which 
identifies  regions  in  phase  space  where  a  transition  can  occur.  It  describes  the  flow  from  each  part  of  phase 
space  to  another,  and  as  the  standard  deviation  of  the  noise  is  increased,  flow  occurs  between  basins.  Using 
the  stochastic  operator,  we  can  construct  a  transition  matrix,  which  yields  the  probability  of  transition 
between  basins  due  to  noise.  By  identifying  the  regions  where  the  solution  has  the  highest  probability  of 
switching  states,  controlling  the  dynamics  may  be  realized. 

2.  Stochastic  Prediction  and  Control 

Predictability  of  seasonally  driven  diseases  that  are  stochastic  is  necessary  for  the  application  of  methods  to 
suppress  future  outbreaks.  Many  vaccine  schemes  are  available  for  equilibrium  diseases,  but  in  the  case  of 
non-equilibrium  outbreaks,  current  methods  may  enhance  outbreaks  or  fail  to  achieve  their  goals.  To 
address  the  problem  of  suppressing  outbreaks  in  stochastic  epidemics,  we  apply  a  new  mathematical  method 
to  a  stochastic  model  to  predict  outbreaks  before  they  occur,  and  then  adapt  a  vaccine  strategy  that  prevents 
the  outbreak  from  occurring.  The  theory  exploits  a  transition  probability  description  from  small  amplitude 
incidence  to  outbreak  dynamics,  and  generates  a  region  of  high  probability  transport  of  the  most  sensitive 
regions  to  stochastic  effects.  Moreover,  it  allows  us  to  monitor  regions  of  stochastic  dynamics  that  have  a 
high  probability  of  preceding  a  large  outbreak,  which  in  turn  leads  to  a  design  of  a  vaccine  control  strategy 
to  suppress  outbreaks.  We  have  developed  a  simple,  but  effective,  general  control  technique  that  takes 
advantage  of  complicated  interactions  of  determinism  and  noise. 


3.  A  noise  induced  bifurcation  observed  experimentally 

The  PI  has  completed  a  successful  collaboration  reproducing  a  3-dimensional  stochastic  bifurcation  both 
experimentally  and  in  a  theoretical  model.  The  physical  model  is  that  of  a  class  B  laser,  which  is  perturbed 
stochastically.  The  effect  of  the  noise  perturbations  on  the  dynamics  is  shown  to  change  the  qualitative 
nature  of  the  dynamics  experimentally  from  a  stochastic  periodic  attractor  to  one  of  chaos-like  behavior.  To 
analyze  the  qualitative  change,  we  apply  the  technique  of  the  stochastic  Frobenius-Perron  operator  to  a 
model  of  the  experimental  system.  The  main  result  is  the  identification  of  a  global  mechanism  to  induce 
chaos-like  behavior  by  adding  stochastic  perturbations  in  a  realistic  model  system  of  an  optics  experiment. 

In  quantifying  the  stochastic  bifurcation,  we  have  computed  a  transition  matrix  describing  the  probability  of 
transport  from  one  region  of  phase  space  to  another,  which  approximates  the  stochastic  Frobenius-Perron 
operator.  This  mechanism  depends  on  both  the  standard  deviation  of  the  noise  and  the  global  topology  of  the 
system.  Our  result  pinpoints  regions  of  stochastic  transport  whereby  topological  deterministic  dynamics 
subjected  to  sufficient  noise  induces  chaos-like  behavior  in  both  theory  and  experiment. 

4.  Master’s  thesis  on  coupled  populations  (in  progress) 

This  thesis,  by  Kristen  Viz,  seeks  to  answer  questions  about  disease  dynamics  and  vaccination  plans.  We  are 
concerned  with  how  a  disease  is  spread  through  and  between  two  populations  that  have  seasonal  contact 
rates,  or  time-varying  parameters.  This  has  applications  to  understanding  the  behavior  of  disease  spread 
between  a  metropolitan  city  and  its  suburbs  or  between  two  neighboring  countries  with  different  vaccination 
rates.  The  importance  of  mathematically  simulating  disease  spread  allows  us  to  look  into  the  future,  possibly 
forecasting  potential  problems  as  parameters  change.  We  will  seek  to  determine  how  the  dynamics  of  a 
larger  population  drives  a  smaller  population,  and  vice  versa.  Also,  we  consider  how  the  sizes  of  these 
populations  effect  this  relation.  We  will  determine  whether  our  model  qualitatively  captures  the  dynamics  of 
the  current  measles  epidemic  in  Cameroon.  If  we  can  show  it  accurately  reflects  the  actual  data,  we  will  be 
able  to  mathematically  manipulate  the  parameters  in  our  model  to  suggest  a  method  to  minimize  the 
numbers  of  measles  cases  in  Cameroon.  This  work  is  in  collaboration  with  Derek  Cummings  and  Donald 
Burke  from  Johns  Hopkins  Bloomberg  School  of  Public  Health  and  Ira  Schwartz  from  the  US  Naval 
Research  Laboratory. 

5.  Dengue  collaboration 

Multi-strain  viruses  are  a  significant  threat  because  of  their  relatively  quick  evolution  and  complex 
spreading  mechanisms.  The  PI  has  been  invited  to  participate  in  an  infectious  disease-modeling  consortium 
headed  by  Dr.  Donald  Burke,  Professor  of  International  Health  and  Epidemiology  at  the  Bloomberg  School 
of  Public  Health  at  Johns  Hopkins  University.  Based  on  data  from  850,000  cases  of  dengue  hemorrhagic 
fever  (DHF)  occurring  in  72  provinces  of  Thailand  for  the  period  1983-1997,  the  Burke  group  has 
discovered  a  spatial-temporal  traveling  wave  in  the  incidence  of  DHF  emanating  from  Bangkok  using 
empirical  mode  decomposition.  Building  on  this  work,  the  PI  has  collaborated  to  develop  and  analyze  a  full 
four-serotype  dengue  model  for  that  region  with  antibody-dependent  enhancement  (ADE).  Using  realistic 
assumptions  and  precise  parameters  derived  from  data,  we  have  analytically  proven  the  existence  of  a  Hopf 
bifurcation  as  we  vary  the  ADE  factor,  and  have  numerically  proven  the  bifurcation  to  chaos.  We  are 
currently  working  to  validate  the  results  of  this  new  model  by  information  provided  by  the  data.  Next,  we 
will  extend  these  results  to  multi-patch  networks,  which  add  spatial  information. 


6.  Undergraduate  projects 

The  PI  has  completed  advising  three  undergraduate  projects.  These  students  investigated  problems 
including:  new  seasonal  models  for  malaria,  the  effects  of  time  varying  vaccination  programs,  and  the 
impact  of  immigrant  subpopulations  with  lower  vaccination  coverage.  The  students  participated  in  a  group 
meeting  once  a  week  to  present  useful  analytic  and  numerical  tools  and  a  summary  of  their  current  results. 
They  presented  their  research  as  posters  at  the  Undergraduate  Sigma  Xi  conference  at  St.  Joseph’s 
University  and  as  talks  at  the  Montclair  State  University  Sigma  Xi  student  research  conference. 

7.  Email  virus  spread  on  a  network 

In  this  time  of  rapidly  growing  computer  speed  and  size,  the  public's  need  for  network  and  Internet 
connections  is  insatiable.  New  viruses  can  spread  unchecked  and  infect  large  numbers  of  computers  very 
quickly.  A  simple  defensive  starting  point  would  be  to  gather  data  that  could  determine  whether  the  virus  is 
in  the  growth  or  die  out  phase,  but  complete  data  sets  for  the  Internet  are  difficult,  if  not  impossible,  to 
obtain.  So,  analysis  must  use  random  samples  and  identify  general  trends.  We  focus  on  viruses  propagated 
by  email.  The  objective  is  to  provide  a  model  that  describes  the  long-term  dynamics  of  email  virus 
propagation  over  a  computer  network.  While  complex,  realistic  models  can  quantitatively  predict  the 
magnitude  of  the  spread,  simpler  ordinary  differential  equations  models  qualitatively  capture  important 
phenomena  like  die  out.  Drawing  on  the  basic  mass-action  models  that  describe  disease  spread  in  biological 
models,  we  use  an  SIS  model  with  parameters  fitted  from  real  data  from  four  viruses.  Therefore,  we  can 
capture  the  current  basic  reproduction  number,  predicting  the  persistence  or  die  out  of  the  virus. 
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Many  mechanical  systems  consist  of  continuum  mechanical  structures,  having  either  linear  or 
nonlinear  elasticity  or  geometry,  coupled  to  nonlinear  oscillators.  In  this  paper,  we  consider  the  class 
of  linear  continua  coupled  to  mechanical  pendula.  In  such  mechanical  systems,  there  often  exist 
several  natural  time  scales  determined  by  the  physics  of  the  problem.  Using  a  time  scale  splitting, 
we  analyze  a  prototypical  structural  -mechanical  system  consisting  of  a  planar  nonlinear  pendulum 
coupled  to  a  flexible  rod  made  of  linear  viscoelastic  material.  In  this  system  both  low-dimensional 
and  high-dimensional  chaos  is  observed.  The  low-dimensional  chaos  appears  in  the  limit  of  small 
coupling  between  the  continua  and  oscillator,  where  the  natural  frequency  of  the  primary  mode  of 
the  rod  is  much  greater  than  the  natural  frequency  of  the  pendulum.  In  this  case,  the  motion  resides 
on  a  slow  manifold.  As  the  coupling  is  increased,  global  motion  moves  off  of  the  slow  manifold  and 
high-dimensional  chaos  is  observed.  We  present  a  numerical  bifurcation  analysis  of  the  resulting 
system  illustrating  the  mechanism  for  the  onset  of  high-dimensional  chaos.  Constrained  invariant 
sets  are  computed  to  reveal  a  process  from  low-dimensional  to  high-dimensional  transitions. 
Applications  will  be  to  both  deterministic  and  stochastic  bifurcations.  Practical  implications  of  the 
bifurcation  from  low-dimensional  to  high-dimensional  chaos  for  detection  of  damage  as  well  as 
global  effects  of  noise  will  also  be  discussed.  ©  2004  American  Institute  of  Physics. 

[DOl:  10.1063/1.1651691] 


Transition  to  chaos  has  been  a  fundamental  problem  in 
nonlinear  dynamics.  The  well  known  routes  to  chaos, 
which  include  the  period-doubling  bifurcation  route,  the 
intermittency  route,  the  quasiperiodic  route,  and  the  cri¬ 
sis  route,  are  for  transition  to  low-dimensional  chaotic 
attractors  with  one  positive  Lyapunov  exponent.  Transi¬ 
tions  to  high-dimensional  chaotic  attractors  with  multiple 
positive  Lyapunov  exponents  have  begun  to  be  addressed. 
Here  we  present  a  class  of  physical  systems  consisting  of 
linear  continuum  mechanical  structures  coupled  to  non¬ 
linear  oscillators.  These  systems  arise  naturally  in  many 
important  engineering  and  defense  applications.  Math¬ 
ematically,  such  a  system  is  typically  described  by  a  set  of 
coupled  partial  and  ordinary  differential  equations, 
which  is  generally  not  amenable  to  analysis.  However,  if 
the  system  exhibits  intrinsically  distinct  time  scales,  ap¬ 
proximations  can  be  made  which  mathematically  reduce 
the  coupled  system  to  a  set  of  ordinary  differential  equa¬ 
tions.  Dynamically,  this  is  equivalent  to  decomposing  the 
motions  into  those  having  slow  and  fast  time  scales,  al¬ 
lowing  for  numerical  and  physical  analyses.  If  the  cou¬ 
pling  between  the  continuum  component  and  the  nonlin- 
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ear  oscillator  is  small,  the  dynamics  can  be  regarded  as 
being  confined  to  a  slow,  approximately  invariant  mani¬ 
fold  exhibiting  low-dimensional  chaos.  Motions  away 
from  the  slow'  manifold  are  typically  fast  in  time,  are 
high-dimensionally  chaotic,  and  they  become  important 
when  the  coupling  is  large  or  when  there  is  noise  present. 
The  system  thus  represents  a  paradigm  for  investigating 
fundamental  phenomena  in  nonlinear  and  stochastic  dy¬ 
namics  such  as  the  transition  to  high-dimensional  chaos 
and  noise-induced  high-dimensional  chaotic  attractors. 
Here  we  shall  demonstrate  that  this  is  so. 


I.  INTRODUCTION 

In  many  mechanical  structures  of  significance,  such  as 
ships,  aircraft,  and  space  vehicles,  there  arise  problems  in 
multi-scale  dynamics  due  to  various  physical  factors.2-5  First 
and  foremost  is  that  many  of  the  structures  we  come  to  de¬ 
pend  upon  are  composed  of  many  sub-structures  covering  a 
wide  range  of  sizes,6  as  in  aircraft  carriers  or  the  space  sta¬ 
tion.  Second,  several  orders  of  magnitude  in  flexibility  may 
be  present,  such  as  a  tether  attached  to  a  satellite/  or  differ¬ 
ent  beam  lengths  in  a  large  truss.8  Such  differences  in  spatial 
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scales  have  recently  led  to  a  number  of  new  multi-scale  nu¬ 
merical  modeling  techniques,  such  as  those  applied  to  the 
finite  element  method.1 

As  a  result  of  such  a  wide  range  of  sizes  and  stiffness 
arising  from  coupled  structures  of  differing  elasticity,  one 
expects  that  there  should  exist  a  range  of  corresponding  dy¬ 
namical  responses  in  frequency.  Therefore,  incorporating  the 
relevant  physics  together  in  a  complex  model  will  generate 
multi-scale  dynamics.  That  is,  there  is  a  set  of  dynamics 
which  may  be  complicated  (such  as  high-dimensional 
chaos),  but  has  definitive  multi-scale  structure.  An  excellent 
example  of  such  a  model  is  the  class  of  driven  coupled  con¬ 
tinuum  mechanical  structures  which  may  be  used  to  explore 
nonlinear  vibrations  in  mechanics. 

In  analyzing  the  dynamics  of  coupled  continuum  me¬ 
chanical  models,  one  must  consider  them  as  spatio-temporal 
systems  which  may  support  a  range  of  dynamical  time 
scales.  Excellent  examples  of  multi-scale  behavior  have  been 
studied  in  the  flexible  spherical  pendulum,9  the  dynamics  of 
a  flexible  beam-oscillator  system,10  and  a  flexible  rod- 
pendulum  system.11  These  examples  demonstrate  the  idea 
that  in  addition  to  a  temporal  splitting  between  fast  and  slow 
time  scales,  there  also  corresponds  a  geometric  splitting. 
Multi-scale  behavior  in  mechanics  is  convenient  since  the 
model  in  many  instances  may  be  decomposed  into  a  global 
singular  perturbation  problem.12  Based  on  a  well-developed 
theory,  one  may  construct  rescaled  systems  for  which  the 
dynamics,  under  suitable  hypotheses,  may  reside  on  an  in¬ 
variant  manifold.13  Physically,  this  might  occur  if  one  struc¬ 
ture  is  almost  perfectly  rigid  (fast  time  scale)  and  attached  to 
a  flexible  structure  (slow  time  scale).  An  excellent  example 
is  the  class  of  ’Mast  slow”  continuum  systems  which  con¬ 
sists  of  “soft”  structures  coupled  to  “stiff”  structures. 

This  class  of  “soft  stiff”  engineering  structures  can 
have  very  complicated  dynamics.  Since  the  stiff  part  of  the 
problem  may  be  considered  an  approximation  to  a  perfectly 
rigid  body,  it  is  reasonable  to  assume  that  part  of  the  com¬ 
plexity  originates  within  the  soft,  flexible  structure.  The  geo¬ 
metric  splitting  yields  invariant  manifolds  which  may  in¬ 
clude  chaos  within  the  soft  structure.14  On  the  other  hand,  for 
critical  parameter  choices,  the  dynamics  may  leave  the  mani¬ 
fold,  and  sample  the  rest  of  the  phase  space,  generating  a 
dimension  changing  bifurcation  which  includes  both  fast  and 
slow  time  scales.14 

Since  multi-scale  engineering  structures  may  have  a  di¬ 
mension  changing  bifurcation,  techniques  for  statistically 
quantifying  the  dynamics  in  space  are  needed.  One  such 
powerful  method  is  that  of  the  method  of  snapshots,  based  on 
the  proper  orthogonal  decomposition  (POD),  or  Karhunen 
Loeve  (KL)  techniques.15’16  First  introduced  to  handle  fluid 
dynamics,  these  methods  have  been  successful  in  quantifying 
the  dynamics  in  fluid  structure  interactions,17  spatio- 
temporal  feedback  control,18  nonstationary  flow  transition 
problems,19  and  aerodynamics  foils.20  They  have  been  used 
to  quantify  a  dimension  change  bifurcation  explicitly  in  a 
soft  stiff  system  operating  near  a  resonant  condition.  In 
quantifying  the  dynamics  of  a  spatio-temporal  system,  the 
KL  technique  is  a  powerful  tool  to  describe  the  modal  struc¬ 
ture,  not  only  analytically  but  computationally  as  well.  In 


contrast,  if  one  wishes  to  compute  dynamical  dimension, 
such  as  Lyapunov  dimension,21  the  linear  variational  equa¬ 
tions  need  to  be  solved  along  the  trajectory,  which  is  prohibi¬ 
tively  expensive  even  for  a  modest  system  of  ordinary  dif¬ 
ferential  equations  (ODEs). 

Although  the  KL  methods  are  useful  for  quantifying  a 
dimension  change  in  dynamics,  they  do  not  necessarily  ex¬ 
plain  the  underlying  cause  of  the  bifurcation.  In  particular, 
when  we  think  of  a  dimension  bifurcation,  we  think  of  a 
change  in  dimension,  abrupt  or  continuous,  as  a  parameter  is 
changed.  Normally,  the  parameter  is  changed  deterministi¬ 
cally,  resulting  in  a  change  in  bifurcation  structure.  This  may 
be  explained  as  sufficiently  sampling  parts  of  an  attractor  off 
of  an  invariant  manifold.  Such  bursting  is  typically  observed 
to  be  chaotic,  which  arises  from  an  underlying  deterministic 
chaotic  saddle.22  However,  another  cause  of  dimension 
changing  chaos  may  be  stochastic.  That  is,  if  sufficient  noise 
is  added  to  a  low-dimensional  attractor,  it  may  generate  a 
high-dimensional  noise  induced  attractor  with  a  new  positive 
Lyapunov  exponent.  Noise  induced  chaos,  when  produced 
along  with  a  dimension  changing  bifurcation,  could  cause 
multi-scale  behavior  in  continuum  mechanics.  Associated 
with  noise  induced  chaos  is  the  idea  of  unstable  dimension 
variability  (UDV). 

Unstable  dimension  variability  is  the  changing  of  the 
number  of  local  unstable  directions  along  a  typical  trajectory. 
Mathematically,  it  can  be  described  in  terms  of  how  a  system 
violates  the  properties  of  hyperbolicity.  This  nonhyperbolic- 
ity  has  been  shown  to  be  fundamental  to  chaotic  dynamics, 
particularly  for  the  problem  of  shadowing  of  numerical  tra¬ 
jectories  in  higher  dimensions.23-30  In  Ref  31,  we  reported 
on  noise  induced  chaos  in  a  preliminary  mechanics  example. 
There,  noise  was  used  to  excite  a  positive  Lyapunov  expo¬ 
nent. 

In  this  paper,  wc  report  on  the  status  of  dimension 
changing  bifurcations  in  both  deterministic  and  stochastic 
mechanical  systems.  As  such,  some  of  the  material  is  neces¬ 
sarily  review.  The  paper  is  laid  out  as  follows:  In  Sec.  II,  the 
full  model  of  a  rod  pendulum  system  is  derived  as  a  nonlin¬ 
ear  coupled  PDE  ODE  system.  A  Galerkin  projection  is 
done  to  put  the  model  in  terms  of  an  infinite  system  of 
ODE's.  and  then  a  finite  dimensional  model  is  extracted  to 
study  the  dynamics.  In  Sec.  Hi,  the  bifurcation  structure  is 
presented  for  the  deterministic  systems  derived  in  Sec.  II. 
Evidence  of  UDV  is  presented  in  the  continuation  diagrams. 
Section  III  C  explores  dimension  change  bifurcations  based 
on  KL  methods,  as  well  as  Lyapunov  spectra.  Section  IV 
explores  the  effects  of  noise  on  the  dynamics. 

II.  DYNAMICAL  CONTINUUM  MECHANICS 

In  general,  the  problems  we  consider  here  model  linear 
continua  coupled  to  nonlinear  oscillators.  That  is,  the  prob¬ 
lem  class  is  that  of  linear  PDF/s  which  are  coupled  to  one  or 
more  nonlinear  oscillators  represented  by  ODE’s.  Such 
linear-nonlinear  coupled  systems  are  ubiquitous  in  many  ap¬ 
plications,  and  are  observed  to  exhibit  nonlinear  vibrations  in 
experiments.10  In  this  section,  we  restrict  ourselves  to  models 
of  linear  elastica  in  one  spatial  dimension.  Such  examples 
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include  cantilevered  beams  and  extensible  rods.  In  general,  if 
we  let  W(^t)  denote  a  measure  of  displacement  as  a  func¬ 
tion  of  space  (£)  and  time  (/),  and  let  xp((j)  be  a  forcing 
function,  then  the  general  equations  of  motion  may  be  rep¬ 
resented  as 

Lpme,t)=xp(€,t), 

d2e  d6  (l) 

jP  +  V  +  G(Wt„)]sme+v-fi  =  0, 

plus  the  appropriate  boundary  conditions.  In  Eq.  (I),  0  de¬ 
notes  the  angular  position  of  an  attached  pendulum  at  a  free 
end  of  the  elastica.  The  term  G(  W  u)  is  the  additional  accel¬ 
eration  imposed  on  the  pendulum  from  the  structure.  Since 
there  is  an  external  driving  body  force  on  the  structure,  the 
function  G(W  tt)  will  also  contain  a  time  varying  source, 
which  will  in  general  depend  on  another  oscillator,  such  as  a 
mechanical  shaker  or  periodic  electric  potential. 

The  differential  operator,  L p,  is  assumed  to  be  linear, 
and  depends  on  a  parameter  which  is  a  measure  of  spectral 
splitting  of  the  relevant  time  scales.  For  a  cantilevered  beam, 
it  has  the  form 

LiiW=Ix2k]W^ii+  W'Mi+lUtiWrtK,  (2) 

while  for  an  extensible  rod  (detailed  below) 

LmW=ii2  W%u  -W&+2  ZrtL  .  (3) 

Normally,  the  external  drive  is  decoupled  from  the  rest  of  the 
structure.  That  is,  it  is  assumed  that  the  external  drive  is 
one-way  coupled  to  the  structure.  In  this  paper,  we  allow  the 
frequency  dependence  of  the  drive  to  depend  weakly  on  the 
dynamics  of  the  structure  itself. 

A.  Full  PDE-ODE  system 

In  formulating  the  dynamics  of  such  a  mutually  coupled 
system,  we  follow  Refs.  13  and  22  in  formulating  in  detail  a 
system  based  on  Eq.  (3).  We  consider  a  specific  mechanical 
system  consisting  of  a  vertically  positioned  viscoelastic  lin¬ 
ear  rod  of  density  pr ,  with  cross-section  Ar  and  length  Z.,. , 
with  a  pendulum  of  mass  Mp  and  arm  length  Lp  coupled  at 
the  bottom  of  the  rod  and  where  the  rod  is  forced  from  the 
top  harmonically  with  frequency  D.  and  magnitude  cx.]*  The 
rod  obeys  the  Kelvin-Voigt  stress-strain  relation'2  and  Er 
and  Cr  denote  the  modulus  of  elasticity  and  the  viscosity 
coefficient.  Cp  is  the  coefficient  of  viscosity  (per  unit  length) 
of  the  pendulum  and  g  is  the  gravitational  constant  of  accel¬ 
eration.  The  pendulum  is  restricted  to  a  plane,  and  rotational 
motion  is  possible.  The  system  is  modeled  by  the  following 
equations: 


FIG.  1.  Rod-pendulum  configuration. 


and  where 

Tp  =  MpLp02  +  Mp(g- xA  -  iiB)  cos(  6) , 

denotes  the  tension  acting  along  the  rigid  arm  of  the  pendu¬ 
lum.  The  variable  u(x,t)  denotes  the  displacement  field  of 
the  uncoupled  rod  with  respect  to  the  undeformed  configu¬ 
ration  at  equilibrium,  relative  to  the  point  A,  while  uH  de¬ 
notes  the  relative  position  of  the  coupling  end  B  of  the  rod 
with  respect  to  point  A .  See  Fig.  1  for  a  schematic  of  the  rod 
and  pendulum  system. 

We  further  suppose  that  the  drive  at  Af  given  by  the 
function  xA(t)  in  Eq.  (4),  is  such  that  it  comes  from  another 
oscillator.  We  suppose  that  the  oscillator  is  weakly  coupled 
to  the  pendulum  through  its  frequency.  Specifically,  we 
model  the  drive  oscillator  by 

Ct>  ,  =  <J>  ,  +  U(  1  +  SP(«(X,0))^>2  -  ®  I  (<J>  I  +  (l>2) 

4>2  =  -n(l+2P(i/O,0))4>i  +  «>2-<l,2(^?+<J)2) 

=  F2(<I>i  (5) 

where  P  is  a  projection  onto  a  Fourier  mode  (see  below), 
and  |2|<1  is  the  coupling  term  that  modulates  the  fre¬ 
quency.  Notice  that  when  2  =  0,  the  solution  of  Eq.  (5)  con¬ 
sists  of  sines  and  cosines  of  frequency  o>  given  the  appropri¬ 
ate  initial  conditions.  In  terms  of  the  solutions  to  Eq.  (5), 
note  that  ^(/)“d>2(7,2). 

Equations  (4)  and  (5)  are  nondimensionalized  by  the  fol¬ 
lowing  variable  rescalings: 


MpLp9+  Mp\g-xA  -  i7/j]sin(  0)  +  CpLp0=  0, 
A,.p,.ii(.x.l)-ArEl.u"(x>t)-A,.Crii"(x,t) 

-A,.p,.(g-xA)  =  0,  (4) 

where  =  d/dt,  and  with  boundary  conditions 

du  diiii 

u(x  =  0,/ )  =  0,  A  ,.Er  —  -  A  rEr  — -  =  Tp  cos(  0) , 

ox  .  ox  ' 
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ip 


J_  £jl 

2  <.op  M  p  ’ 


1  7T2  Cr 


where 


77(2///-  1) 

~  Lr 


m 


are  the  natural  frequency  of  the  uncoupled  pendulum  and  the 
spectrum  of  natural  frequencies  of  the  uncoupled  flexible 
rod,  respectively,  while  and  £,.  denote  their  damping  fac¬ 
tors. 

Using  the  parameter  rescalings,  and  setting  time  deriva¬ 
tives  equal  to  zero,  the  stable  and  unstable  static  equilibrium 
configurations  of  the  coupled  rod  and  pendulum  system  are 
given  by  ( 0C,U )  and  (ftv  +  ,(7),  where 


0C=O,  0S  >  =  ±  77, 

„  fl2TT  2 

£/=  ^[2(1 +£)£-**]. 


The  normalized  equations  are  thus 

0+  [  1  -  VB(  t)-XJ  t) lsin(  0)  +  2 6=  0, 

At2 7T2  V( r)  -  r ( t.  r)  -  8  K"(  g,  t)  =  ~V  tt2V,( t), 

F(£=(),t)  =  0,  \,t)=~h2Ptt2[1  -Tcos{0)}, 

(6) 

where 


F(  f ,  r)  -  U( &  r)  -  £),  0^  ^  1 ,  -  oo  <  r<  +  oo, 


Finally,  consider  the  finite  set  of  ordinary  differential 
equations  obtained  from  Eq.  (7)  by  truncating  to  the  first  N 
rod  modes  and  applying  the  additional  rescalings  {'V  { ,XV2} 
=  {0,0}  and  {n2(i2mZ2m-i  =  obtaining 


^7  = 


1-2  (-1 1''  1/;V(^.Z)-«^4 


J= ' 


sin(AJ/'  i ) 


+  2£/I'2, 

^4  =  F2(^3,^4,S.n), 


(8) 


^2m  > 

~  f  > M  “1,2,...  ,N, 


where 


1 


--Z2;/)„t  +  2^Z2;w 


L 

-  ( -  1 )"' :  1 2  p[ 'Vl  cos(  ¥ , )  -  sin2 OP , )] 

[4  M 


77 


■  +  (-l)m+12yBcos2(t|) 


and  note  that  we  redefine  =  B/ dr  and  '  -  for  the  re¬ 
mainder  of  the  paper. 


B.  Projection  onto  a  finite  model 


In  carrying  out  our  analysis,  we  will  consider  a  reduction 
of  the  ODE  PDF  system  in  Eq.  (6).  This  reduction  is  ob¬ 
tained  by  performing  a  modal  expansion  of  the  rod  equation, 
where  the  displacement  V  is  expanded  as  V(  £,  r) 
=  2*^,7 ]m(r) This  results  in  an  infinite  system  of 
coupled  oscillators, 


6 =  ■ 


i  +  2  (-i y'vrx.,(r) 


sin(  0)-2£p0, 


Vm  Vm 

Lm(0)vi=-—T-T  +  Hr - - 

4  rj-  T)-m  fi/J.-, 


-  (  -  I )"'  ’  1 2 /3[  0-  cos(  0)  -  sin2(  0)} 


4 

- +  (- 

77 


)"'i  1 2  /?  cos2(  0) 


XA  (r). 


(7) 


equivalent  to  Eq.  (4).  where  0)  is  the  infinite  linear  op- 
era  tor 


2  r-5,,,,  +  1-  I  >2 fi  cos,2 ($)}. 

j  ~  1 

See  Ref.  22  for  the  details  of  this  transformation. 


^3,  ^4  are  drive  variables  as  defined  in  Eq.  (5),  and 
L~\r(9)  is  the  inverse  of  the  NXN  truncation  of  operator 
Lm{0).  Fj  and  F2  are  given  by  the  right-hand  sides  of  Eq. 
(5).  Note  that  Eq.  (8)  is  an  autonomous  system,  and  the  cy¬ 
clic  variables,  AT3  and  T'4  are  introduced  to  account  for  the 
periodic  forcing,  which  has  period  12  when  the  coupling  pa¬ 
rameter  X  =  0.  For  this  paper,  unless  otherwise  noted,  we 
consider  the  truncated  system  obtained  by  taking  N-  1.  We 
also  considered  the  system  in  Eq.  (8)  with  N=2  and  N 
=  10  and  found  qualitatively  similar  dynamics.  Notice  that 
the  terms  Z2m-\  correspond  to  the  rod  displacement  ampli¬ 
tudes,  while  the  even  indexed  terms  Z2m  are  the  rod  velocity 
mode  amplitudes.  The  function  fN{yV,Z)  is  similar  to  the  one 
defined  in  Ref.  13  in  the  case  where  £  =  0,  and  the  derivation 
may  be  found  there. 

The  primary  parameter  governing  the  coupling  between 
the  rod  and  pendulum  is  the  ratio  of  the  natural  frequency  of 
the  pendulum  to  the  frequency  of  the  first  rod  mode,  ii 
=  o)p!(x)x .  In  the  limit  w j-**,  the  rod  is  perfectly  rigid,  fx 
— >0,  and  the  system  reduces  to  a  forced  and  damped  pendu¬ 
lum.  For  0<ix<^  1  sufficiently  small,  global  singular  pertur¬ 
bation  theory7  predicts  that  system  motion  is  constrained  to  a 
slow  manifold,  and  the  (fast)  linear  rod-modes  are  slaved  to 
the  slow  pendulum  motion.12  For  nonzero  a  (the  amplitude 
of  the  periodic  forcing)  the  slow  manifold  is  a  nonstationary 
(periodically  oscillating)  two-dimensional  surface. 

For  our  study,  we  set  the  number  of  modes  in  tire  struc¬ 
ture  to  be  unity,  and  consider  the  following  reduced  system: 
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Xp2  =  ( 7ry2(  - 2 Cp-4PCP  cos2( ¥ , )  +  p  sin(2'If ,  )'I'2) 

—  sin('P  j )  (4  a'V  A  +  4  a/3  -  cos2  ( |  )'P4 
+  ir(  1  +2p+Zl+2CrZ2+a'VA)))/(TrS), 

^j=^3(i  -(^f+^J))+nn  +  «r^2)^4, 

(9) 

^4  =  ’Jr4(  1  -  ('?!  +  ^4»  -  fi  (  1  +  <r^2>^3  , 

Z}  =Z2  //LI, 

Z,  =  -  (2/?tt-  2/?tt  cos('I' ,  )V\  +  irZ,  +  2  tt£,Z2 
+  4a1Ir4  +  2/8'7rcos2('Ir1)(—  1  +«'Ir4))/(/4.'!r5), 

where  8~  1  +2/? cosfT'j).  In  Eq.  (9),  XV  x  and  XV2  are  the 
pendulum  position  and  velocity,  XV  $  and  'I'4  are  the  drive 
oscillator  variables,  and  Zx  and  Z2  are  the  rod  mode  position 
and  velocity.  Notice  that  the  singular  perturbation  parameter 
denotes  the  rod  variables  to  have  a  fast  time  scale  compared 
to  the  pendulum  and  drive.  Also,  the  frequency  of  the  drive 
oscillator  is  a  function  of  the  pendulum  momentum,  which  is 
a  feedback  term.  The  system  in  Eq.  (9)  is.  therefore,  fully 
coupled,  which  is  a  generalization  of  the  more  ideal  case  of 
having  a  perfectly  isolated  drive. 


111.  BIFURCATION  STRUCTURE  OF  THE 
DETERMINISTIC  SYSTEM 

We  now  consider  the  deterministic  one  rod  mode  model 
obtained  in  Sec.  II.  It  is  useful  to  first  consider  the  determin¬ 
istic  model,  without  any  added  noise.  The  underlying  dy¬ 
namical  structures  of  the  deterministic  system  will  determine 
in  what  way  additive  noise  manifests  itself  in  the  dynamics. 
While  Eq.  (9)  is  much  less  complex  than  the  original  PDE,  it 
still  exhibits  a  wide  variety  of  complicated  behaviors.  In  par¬ 
ticular,  when  the  amplitude  a  of  the  forcing  is  sufficiently 
large,  solutions  of  Eq.  (9)  are  chaotic,  and  such  solutions 
with  both  one  and  two  positive  Lyapunov  exponents  have 
been  observed.  In  addition  to  the  forcing  amplitude  a,  the 
behavior  of  solutions  of  Eq.  (9)  is  dramatically  affected  by 
the  value  of  the  coupling  parameter  (i. 

Since  Eq.  (9)  is  singularly  perturbed  for  0 <  /x 1 ,  we 
will  obtain  a  description  of  the  slow  dynamics  by  closely 
following  the  geometric  approach  adopted  in  Ref. 
22.  The  slow  manifold  approximation,  {Zi,Z2} 
=  7/m(vP j  ,xT2  AF3,'vR4)  is  obtained  for  Eq.  (9)  (with  cr=0), 
using  the  method  of  Ref.  12,  The  slow  manifold  (given  by 
the  graph  of  H^)  is  a  submanifold  in  phase  space  on  which 
the  slow  dynamics  reside,  and  which  relates  the  rod  motion 
to  the  pendulum  motion  and  periodic  forcing.  When  jjl  is 
sufficiently  small,  the  dynamics  of  Eq.  (9)  can  then  be  ap¬ 
proximated  by  the  reduced  system: 


^,  =  ^2. 

^2  =  ( - 2 ip- ~ 4 PCP  cos2 OP , )  +  p  sin( 24',)'T'2) 

—  sin(  | )  ( 4  4  +  4  7r  cos2  ( 'P  |)1P4 

+  Tr(\+2p+H/J.(-)  +  2£rHli(-)  +  a'V4)))/(ir8), 

,  (10) 
'p3='p3(i-('p5+'p4))+n(i+o-'p2)'p4, 

'P4  =  'P4(  1  -  (^1+  ¥4)) ~fi(  1  +  o-Mr2OP3 . 

It  is  important  to  note  that  the  above  approximation  is 
valid  only  for  /x  sufficiently  small.  In  particular,  it  has  been 
noticed  that  when  fi  is  increased,  the  dynamics  of  the  full 
one  rod-mode  system  given  by  Eq.  (9)  no  longer  remains  on 
the  slow7  manifold,  but  exhibits  a  bursting  characteristic.  This 
bursting  is  not  directly  observed  in  the  system  variables 
themselves,  but  rather  when  observing  the  variable  defined 
by  A  =  N/(Z, -  ))2  +  (Z2  -//„,2( •  W,  where 

denotes  the  / th  component  of  H ^ .  The  variable  A  measures 
the  distance  of  the  solution  of  the  rod  components  of  the  full 
problem  to  the  slow  manifold.  For  /x  sufficiently  small,  only 
small  bursts  A  are  observed.  Elowever,  as  fx  is  in¬ 

creased,  the  amplitude  of  bursts  quickly  increases. 

Examining  the  variable  Aj^Zj -H^\(  • )  gives  the  dif¬ 
ference  between  the  actual  rod  displacement  and  the  slaved 
rod  displacement  as  calculated  using  the  slow7  manifold  ap¬ 
proximation.  As  /x  increases  so  that  A  is  0(  1 )  in  amplitude, 
A 1  has  the  appearance  of  a  relaxation  oscillation.  That  is,  a 
large  excursion  can  be  observed  away  from  the  slow7  mani¬ 
fold  approximation  with  succeeding  bursts  decaying  in  am¬ 
plitude  toward  the  slow  manifold,  and  another  large  burst 
may  occur  before  the  solution  has  reached  the  slow7  mani¬ 
fold. 

For  small  /x,  solutions  of  Eq.  (10)  agree  well  with  solu¬ 
tions  of  the  full  system  modeled  by  Eq.  (9).  This  ceases  to  be 
the  case  as  /x  is  increased,  and  nontrivial  fast  dynamics  de¬ 
velop.  We  will  consider  two  particular  choices  of  the  cou¬ 
pling  parameter,  to  illustrate  the  nontrivial  fast  dynamics 
which  develop.  In  the  next  subsection  we  examine  the  dy¬ 
namics  of  Eq.  (9)  when  the  coupling  is  relatively  small:  fx 
=  0.086  875.  For  this  value  of  /x,  system  motion  is  no  longer 
confined  to  the  slow7  manifold.  Instead  bursting  off  the  slow 
manifold  is  observed.  We  will  then  examine  the  dynamics  of 
Eq.  (9)  for  the  relatively  large  value  of  jx~  0.5025  in  the 
following  subsection.  Near  this  value  of  /x  there  is  an  internal 
2:1  resonance,  and  the  observed  bursting  is  much  greater  in 
amplitude.  For  the  simulations  presented  in  the  remainder  of 
this  section,  we  set  (3=  1,  £,=  £^  =  0.01,  0.0001,  and  f! 

=  1.9527  unless  otherwise  stated. 

A.  Lyapunov  exponents  and  attractor  bifurcation 
structure:  The  singularly  perturbed  case 

Recall  that  as  fi-> 0,  there  ceases  to  be  any  nontrivial 
rod  motion  whatsoever,  and  Eq.  (9)  reduces  to  a  forced  and 
damped  pendulum.  Increasing  /x  amounts  to  increasing  the 
flexibility  of  the  rod,  and  nontrivial  rod  motions  independent 
of  the  pendulum  motion  are  observed.  For  the  simulations  of 
this  subsection,  w7e  set  ^  =  0.086  875,  and  vary  cr,  while  the 
other  parameters  are  as  stated  in  the  introduction  of  this  sec- 
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H=0.086875,  <5=0.0001 


FIG.  2.  The  attractor  bifurcation  diagram  for  /x~ 0.086  875  and  a 
=  0.0001,  showing  the  rod  displacement  A,=Zt  -  \  ,^2»^3  >^4)  as 

a  function  of  the  forcing  amplitude  a.  There  is  a  stable  periodic  orbit  for 
1.30. 


tion  above.  We  choose  this  value  for  fi  because  it  is  large 
enough  that  there  is  significant  motion  on  fast  time  scales, 
but  small  enough  that  our  slow  manifold  approximation  will 
be  reasonably  accurate. 

We  consider  first  the  attractor  bifurcation  structure  of  Eq. 
(9).  In  particular,  we  will  use  the  bursting  variable  Aj  as  a 
function  of  a .  We  examine  the  Poincare  map  defined  by 
strobing  the  flow  whenever  3  =  0  and  XV4  =  1.  For  a  suffi¬ 
ciently  small  (A1=0)  the  pendulum  is  stationary  and  only 
periodic  motions  on  the  fast  manifold  are  observed.  As  a 
increases  through  0.17.  the  pendulum  transitions  to  periodic 
motion,  and  the  resulting  periodic  orbit  slaves  the  system  to 
the  slow  manifold.  For  a**  1 .29,  there  is  the  sudden  onset  of 
chaos.  See  Fig.  2. 

To  further  elucidate  the  bifurcation  structure  of  the  one 
mode  model,  we  utilized  the  bifurcation  continuation  soft¬ 
ware  AUTO  97, 33  using  a  as  our  continuation  parameter.  We 
started  the  continuation  calculation  using  a  periodic  orbit  we 
calculated  numerically  for  a- 0.01.  For  this  value  of  a,  we 
found  the  stable  periodic  orbit  to  consist  of  a  motionless 
pendulum  and  an  oscillating  rod  mode.  Physically,  this  cor¬ 
responds  to  the  pendulum  hanging  straight  down,  with  the 
only  motion  consisting  of  deformations  of  the  rod  with  the 
same  period  as  the  periodic  forcing.  There  is  a  secondary 
branch  which  is  born  in  a  saddle-node  bifurcation  at  a 
0.044.  The  saddle  branch  meets  the  first  branch  of  nodes  at 
a  period  doubling  bifurcation,  while  the  upper  branch  of 
nodes  ends  at  a  torus  bifurcation  at  1.29.  This  agrees 
closely  with  what  is  observed  in  the  attractor  bifurcation  dia¬ 
gram.  See  Fig.  3. 

We  additionally  computed  the  Lyapunov  exponents  of 
the  one  rod  mode  model,  Eq.  (9),  shown  in  Fig.  4.  The  most 
notable  feature  of  the  Lyapunov  spectrum  for  fi  =  0.086  875 


FIG.  3.  The  continuation  diagram..  The  lower  branch  consists  of  a  solution 
in  which  the  pendulum  does  not  swing  and  only  rod  motion  is  present,  while 
the  upper  branch  corresponds  to  nontrivial  periodic  motion  of  both  the  pen¬ 
dulum  and  the  rod.  Solid  dots  denote  a  stable  periodic  orbit  while  open  dots 
denote  an  unstable  periodic  orbit.  The  upper  stable  branch  exists  for  a 
e  (0.044,1.29).  As  a  increases  through  1.29,  stability  of  the  periodic  orbit  is 
lost,  closely  corresponding  with  the  behavior  seen  in  Fig.  2. 

is  that  for  .75  chaotic  orbits  have  one  positive  Lyapunov 
exponent,  while  for  a S:  1.75,  chaotic  orbits  have  two  posi¬ 
tive  Lyapunov  exponents,  and  the  transition  from  one  to  two 
positive  Lyapunov  exponents  is  smooth.  This  implies  that  as 
the  forcing  amplitude  a  is  increased,  the  chaotic  attractor 
increases  in  dimension. 

Finally,  we  used  the  constrained  invariant  manifold 
(CIM)  method  exposited  in  Ref.  22  to  compute  the  approxi¬ 
mation  of  the  stable  manifold  of  the  chaotic  saddle  con¬ 
strained  to  the  slow  manifold.  Briefly,  the  method  works  by 
finding  those  initial  conditions  on  the  slow  manifold  which 

0  .  1 

remain  within  €+  of  the  slow  manifold  to  time  T  under 
evolution  of  Eq.  (9).  Figure  5  shows  the  projection  of  the 
approximated  manifold  onto  the  pendulum  variables  {  and 
'PV  This  set  gives  an  indication  as  to  the  structure  of  the 


FIG.  4.  The  Lyapunov  spectrum  for  0.086  875.  The  most  negative 
Lyapunov  exponent  is  not  shown. 


Downloaded  02  Jul  2004  to  137. 350.140  107  Redistribution  subject  to  AIR  license  or  copyright,  see  http://chnos.nl;  org/chnos/copyright.jsp 


Chaos,  Vol.  14,  No.  2,  2004 


*! 


FIG.  5.  The  stable  manifold  of  the  chaotic  saddle,  constrained  to  the  slow 
manifold,  calculated  using  the  CIM  method  and  projected  onto  the  pendu¬ 
lum  subspace  The  parameters  used  were  /a=0.086  875  and  a 

=  1.2,  while  /?,  cr,  £r ,  and  £p  are  as  noted  in  the  text. 

invariants  constrained  to  the  slow  manifold.  The  CIM 
method  parameters  used  were  T+  =  8.0  and  €+  =0.1. 

B.  Lyapunov  exponents  and  attractor  bifurcation 
structure:  The  near  resonance  case 

We  next  increase  jx  to  0.5025.  For  this  value  of  the  cou¬ 
pling,  there  is  an  internal  2:1  resonance  in  the  system. 

As  in  the  previous  subsection,  we  consider  the  attractor 
bifurcation  structure  of  Eq.  (9),  measuring  the  bursting  vari¬ 
able  A  j  as  a  function  of  a .  We  again  use  the  Poincare  map 
given  by  strobing  the  flow  whenever  ^3  =  0  and  ^4=I. 
Since  ju  is  relatively  large,  we  do  not  expect  the  slow  mani¬ 
fold  approximation  H jJL(-)  to  be  very  accurate.  Indeed  for  a 
small,  we  observe  that  A,  is  no  longer  zero,  but  the  slow 
manifold  approximation  indicates  the  solution  lies  an  0(1) 
distance  from  the  slow  manifold.  However,  we  still  expect 
that  the  observed  stable  periodic  motions  are  slow.  For  a 
^0.88,  there  is  the  sudden  onset  of  chaos,  and  the  solutions 
of  Eq.  (9)  move  far  from  the  slow  manifold  approximation, 
as  shown  in  Fig.  6. 

Using  a  as  our  continuation  parameter,  we  started  the 
continuation  calculation  using  a  periodic  orbit  we  calculated 
numerically  for  a~  0.01  (see  Fig.  7).  For  this  value  of  a ,  we 
again  find  the  stable  periodic  orbit  to  consist  of  a  motionless 
pendulum  and  an  oscillating  rod  mode.  There  is  a  secondary 
branch  which  is  bom  in  a  saddle-node  bifurcation  at  a 
^0.044.  The  saddle  branch  meets  the  first  branch  of  nodes  at 
a  period  doubling  bifurcation,  while  the  upper  branch  of 
nodes  ends  at  a  torus  bifurcation  at  <*^0.27.  This  agrees 
closely  with  what  is  observed  in  the  attractor  bifurcation  dia¬ 
gram,  as  shown  in  Fig.  6. 

Finally,  we  examine  the  Lyapunov  spectrum  of  Eq.  (9) 
near  resonance  shown  in  Fig.  8.  The  most  notable  feature  in 
this  case,  is  the  sudden  transition  from  a  stable  periodic  orbit 
to  hyperchaos.  In  fact,  in  numerical  studies  on  a  mesh  of  fx 
values  of  width  0.021875  from  yu  =  0.086875  to  /x 
=  0.5025,  we  did  not  observe  a  smooth  transition  from  one 
to  two  positive  Lyapunov  exponents  for  any  //,>0.086 875. 
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jt=0.5025,  g=0,0001 


FIG.  6.  The  attractor  bifurcation  diagram  for  fi  =  0.5025  and  a~  0.0001, 
showing  the  rod  displacement  A]  as  a  function  of  the  forcing  amplitude  a. 
A  stable  periodic  orbit  is  observed  for  a:SO,88. 


Rather,  we  observed  an  apparent  discontinous  change  in  the 
distribution  of  Lyapunov  exponents,  similar  to  the  behavior 
observed  in  Ref.  34  where  the  sudden  transition  to  hypercha- 
otic  behavior  is  observed  in  the  same  model  with  (7=0,  that 
is,  without  feedback  to  the  drive.  The  transition  to  hyper¬ 
chaos  occurs  near  a^0.88,  which  agrees  well  with  the  ter¬ 
mination  of  the  lowest  stable  branch  of  nodes  (see  Fig.  7). 

We  again  apply  the  CIM  method  to  Eq.  (9),  this  time  for 
/a  =  0.5025,  as  shown  in  Fig.  9.  Since  the  system  is  in  reso¬ 
nance  and  we  are  using  the  slow  manifold  approximation  at 
the  edge  of  its  applicability,  we  set  the  threshold  £+  =  150, 


FIG.  7.  The  continuation  diagram.  The  lower  branch  consists  of  a  solution 
in  which  the  pendulum  docs  not  swing  and  only  rod  motion  is  present  while 
the  upper  branches  corresponds  to  nontrivial  periodic  motion  of  both  the 
pendulum  and  the  rod.  Solid  dots  denote  a  stable  periodic  orbit  while  open 
dots  denote  an  unstable  periodic  orbit.  Note  that  the  lowest  stable  branch 
pictured  ends  in  a  saddle-node  bifurcation  (label  4)  at  a^O.88,  correspond¬ 
ing  to  the  point  in  the  attractor  bifurcation  diagram  (Fig.  6)  where  a  chaotic 
solution  is  first  seen. 
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HU.  8.  The  Lyapunov  spectrum  for  /i- 0.5025.  The  most  negative 
Lyapunov  exponent  is  not  shown. 


while  keeping  =  8.  The  large  value  for  £+  may  be  justi¬ 
fied  by  the  fact  that  the  burst  extrema  observed  in  A,  are 
now  much  larger  in  amplitude.  There  is  clearly  an  interesting 
level  of  structure  in  the  resulting  set,  though  just  how  this  set 
should  be  interpreted  remains  an  open  question. 

C.  Dimension  changing  bifurcations 

In  order  to  examine  the  fidelity  of  the  one  rod  mode 
model,  we  examine  the  same  model,  but  this  time  truncating 
the  rod  expansion  at  16  modes.  We  find  that  while  the  main 
attractor  bifurcations  occur  for  smaller  values  of  the  forcing 
amplitude  a,  the  bifurcation  structure  has  a  similar  qualita¬ 
tive  appearance  to  the  one  mode  resonant  case.  Due  to  the 
expense  of  calculating  the  slow  manifold  approximation,  m 
Fig.  10  we  plot  the  first  position  rod  mode  against  the  forcing 


2[ - T - 1 - T - 1 - T - 1  I  r 


FIG.  9.  The  stable  manifold  approximation  of  the  chaotic  saddle,  con¬ 
strained  to  the  slow  manifold,  and  projected  onto  the  pendulum  subspace 
(xp,  ,\p2),  calculated  using  the  C1M  method.  The  parameters  used  were  /x 
=  0.5025  and  n-  =  0.7,  while  ft  <r,  and  are  as  noted  in  the  text. 
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H=0.5025,  a=0.0001 ,  N=16 


FIG.  10.  The  attractor  bifurcation  diagram  for  fi— 0.5025  and  <7  0.0001 
for  the  16  rod  mode  model,  showing  the  rod  displacement  as  a  function  of 
the  forcing  amplitude  a.  A  stable  periodic  orbit  is  observed  for  aS0.32. 


amplitude  a.  The  behavior  of  the  solutions  is  similar  for 
what  we  observe  in  the  one  mode  resonant  case.  There  is  a 
periodic  orbit  that  abruptly  transitions  to  a  chaotic  solution 

near  at  =  0.32. 

We  again  ran  AUTO  on  the  16  mode  model  as  shown  in 
Fig.  1 1 .  The  saddle  structure  is  somewhat  different  from  that 
of  the  one  mode  resonant  case  shown  in  Fig.  7,  but  some  of 
the  essential  features  are  still  observed.  In  particular,  the 
lower-most  stable  branch  terminates  in  a  saddle-node  bifur¬ 
cation  near  a  =  0.3,  close  to  the  value  for  a  at  which  we 
observe  the  first  chaotic  motion,  as  seen  in  Fig.  1 0. 


FIG.  1 1  The  continuation  diagram  computed  using  AUTO.  The  lower  branch 
consists  of  a  solution  in  which  the  pendulum  does  not  swing  and  only  rod 
motion  is  present  while  the  upper  branches  corresponds  to  nontrivial  peri¬ 
odic  motion  of  both  the  pendulum  and  the  rod.  Solid  dots  denote  a  stable 
periodic  orbit  while  open  dots  denote  an  unstable  periodic  orbit.  Note  that 
the  lowest  stable  branch  pictured  ends  in  a  saddle-node  bifurcation  (label  4) 
at  »~0.3,  corresponding  to  the  point  in  the  attractor  bifurcation  diagram 
(Fig.  10)  where  a  chaotic  solution  is  first  seen. 
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FIG.  12.  The  number  of  KL  modes  required  to  capture  99%  of  the  system 
energy,  as  a  function  of  forcing  amplitude  a. 


It  is  expensive  to  calculate  the  Lyapunov  exponents  for¬ 
th  is  system,  since  for  a  system  of  N  equations,  one  is  re¬ 
quired  to  integrate  the  N  differential  equations,  plus  N2 
variational  equations  in  order  to  approximate  the  Lyapunov 
exponents.  For  the  16  mode  truncation,  we  instead  compute 
the  KL  dimension  of  the  attractor  as  shown  in  Fig.  12. 
Briefly,  the  KL  decomposition  gives  the  optimal  way  to  com¬ 
pute  an  orthogonal  linear  expansion,  in  an  energy  sense. 
Thus,  the  first  mode  in  a  KL  expansion  will  contain  the 
maximal  energy  possible  for  a  linear  mode,  the  second  will 
contain  the  second-most  energy  possible,  and  so  on.  Then,  a 
definition  for  KL  dimension  is  the  minimal  number  of  KL 
modes  required  to  satisfy  some  (large)  energy  threshold. 
Briefly,  we  describe  how  to  compute  the  KL  dimension,  and 
the  reader  should  see  Ref.  13  for  details. 

We  consider  the  system  given  by  Eq.  (8),  and  define  the 
field  to  be  the  vector  of  continuous  functions  of  time  defined 
by 

{/(/)  =  ['!' 1  2  3  ^ 4  \  l  yZhvVO)- 

(ID 

Computing  the  KL  modes  is  based  on  a  method  which  maxi¬ 
mizes  the  variance  and  minimizes  the  covariance.  For  two 
time  sampling  (tm  we  define  the  nm  entry  of  the  corre¬ 
lation  matrix  as 

C,"”  =  ^jUT(‘n)U(tm).  (12) 

M  is  the  number  of  time  snapshots,  and  the  indices  n,m 

-  1,2,..., A/. 

The  “energy”  of  the  system  is  quantified  by  measuring 
the  number  of  active  modes,  which  is  done  by  gleaning  in¬ 
formation  from  the  spectrum  of  the  correlation  matrix  C,  i.e., 
by  solving  the  eigenvalue  problem 

CAk=\kAk.  (13) 


The  sum  of  the  Xk  is  the  total  energy  of  the  system,  and  the 
sum  is  used  to  normalize  the  spectrum.  That  way,  we  can 
choose  a  threshold,  say  99%,  and  pick  those  modes  that  are 
in  the  sum.  We  can  then  define  a  KL  dimension  of  the  dy¬ 
namics  that  consists  of  those  modes  that  make  up  99%  of  the 
energy. 

For  a  series  of  values  of  a,  we  computed  the  number  of 
KL  modes  required  to  capture  99%  of  the  system  energy.  We 
find  that  as  a  increases  into  the  chaotic  regime,  the  number 
of  modes  needed  to  capture  the  system  energy  increases  dra¬ 
matically.  This  implies  that  the  underlying  chaotic  solutions 
are  high  dimensional. 

IV.  NOISE  INDUCED  CHAOS 

In  the  previous  sections,  deterministic  bifurcations  were 
considered  in  which  both  low  and  high  dimensional  dynam¬ 
ics  were  observed.  However,  in  many  real  engineering  sys¬ 
tems  of  interest,  noise  plays  a  definitive  role  with  far  reach¬ 
ing  consequences,  as  seen  in  Refs.  35-37.  In  fact,  it  is 
possible  to  generate  chaos  using  additive  noise  in  mechani¬ 
cally  driven  systems,  as  seen  in  the  example  of  driven  sto¬ 
chastic  mechanics  presented  in  Ref.  31.  In  this  section,  we 
wish  to  take  the  liberty  to  quantify  how  noise  induced  chaos 
is  related  to  a  novel  mathematical  quantity  related  to  the 
unstable  dimensionality  of  the  system.  Once  the  dynamical 
systems  are  sufficiently  high  dimensional,  it  will  be  seen 
how  noise  interacts  with  unstable  spaces  to  produce  positive 
Lyapunov  exponents.  Since  flexible  continuum  mechanical 
systems  produce  high  dimensional  dynamics,  the  role  of 
noise  will  be  seen  to  be  play  a  prominent  role  in  bifurcation 
theory. 

A.  Unstable  dimension  variability  associated  with 
noise-induced  chaos  and  scaling  law  of  Lyapunov 
exponents 

1.  Noise-induced  unstable  dimension  variability 

An  interesting  phenomenon  associated  with  noise- 
induced  chaos  is  that  unstable  dimension  variability  arises  as 
soon  as  the  attractor  becomes  chaotic.  Unstable  dimension 
variability  means  that,  along  a  typical  trajectory,  the  number 
of  local  unstable  directions  can  change.  This  is  the  type  of 
nonhyperbolicity  that  has  been  shown  to  be  fundamental  to 
chaotic  dynamics,  particularly  for  the  problem  of  shadowing 
of  numerical  trajectories  in  high  dimensions.23' 30  Math¬ 
ematically,  unstable  dimension  variability  can  be  described 
in  tenns  of  the  notion  of  hyperbolicity  (or  nonhyperbolicity). 

Consider  a  chaotic  set  from  an  /^-dimensional  map.  The 
set  is  hyperbolic  if  the  following  three  conditions  are  met:  lS 
(1)  At  each  point  in  the  set  the  tangent  space  can  be  split  into 
an  expanding  subspace  and  a  contracting  subspace.  Dis¬ 
tances  in  the  expanding  (contracting)  subspace  grow  (shrink) 
exponentially  in  time;  (2)  the  angle  between  the  stable  and 
the  unstable  subspaces  is  bounded  away  from  zero;  (3)  the 
expanding  subspace  evolves  into  the  expanding  one  along  a 
typical  trajectory  and  the  same  is  true  for  the  contracting 
subspace.  Violation  of  condition  (2)  leads  to  nonhyperbolic¬ 
ity  with  tangencies,  which  occurs  commonly  in  low- 
dimensional  chaotic  systems  with  only  one  unstable  direc¬ 
tion.  Nonhyperbolicity  with  unstable  dimension  variability  is 
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caused  by  the  violation  of  condition  (3),  which  occurs  in 
systems  with  more  than  one  unstable  direction,  i.e..  high- 
dimensional  chaotic  systems.  In  high  dimensions,  commonly 
there  are  systems  that  violate  both  conditions  (2)  and  (3). 

We  can  argue  that  when  noise  induces  a  chaotic  attractor, 
unstable  dimension  variability  arises  immediately.  Consider 
the  common  situation  where  there  are  two  coexisting  dy¬ 
namical  invariant  sets  with  distinct  unstable  dimensions.  For 
instance,  in  the  simplest  case  of  a  one-dimensional  map.  in  a 
periodic  window  an  attracting  periodic  orbit  with  zero  un¬ 
stable  dimensions  coexists  with  a  chaotic  saddle  with  un¬ 
stable  dimension  one.  Besides  periodic  windows,  another 
situation  is  where  there  is  a  periodic  attractor  and  several 
isolated  saddle  periodic  orbits.  The  stable  and  unstable  mani¬ 
folds  of  these  orbits  are  close  to  each  other  and  are  about  to 
form  homoclinic  or  heteroclinic  intersections.  The  presence 
of  noise  can  materialize  the  intersections,  creating  a  chaotic 
set,  the  so-called  stochastic  chaotic  saddle . 39-40 

For  any  periodic  point  on  the  attractor,  under  additive 
noise  of  amplitude  D  a  trajectory  can  be  found  in  a  ball  of 
radius  D.  If  D  is  small  so  that  the  ball  does  not  intersect  the 
stable  manifold  of  the  chaotic  saddle,  the  final  attractor  of  the 
system  will  simply  be  a  fattened  version  of  the  original  pe¬ 
riodic  attractor.  This  is  so  because  a  random  initial  condition 
leads  to  a  trajectory  that  is  confined  in  the  vicinity  of  the 
periodic  attractor,  although  there  can  be  transient  chaos  ini¬ 
tially,  in  the  sense  that  the  trajectory  may  move  toward  the 
chaotic  saddle  along  its  stable  manifold,  wander  near  the 
saddle  for  a  finite  amount  of  time,  and  leave  it  along  its 
unstable  manifold.  Assume  that  for  Z)  =  Dr,  the  noisy  ball 
begins  to  intersect  the  stable  manifold  of  the  chaotic  saddle. 
For  D>DC ,  there  is  a  nonzero  probability  that  a  trajectory  in 
the  vicinity  of  the  original  periodic  attractor  is  kicked  out  of 
the  noisy  ball  and  moves  toward  the  chaotic  saddle  along  its 
stable  manifold.  Due  to  the  nonattracting  nature  of  the  cha¬ 
otic  saddle,  the  trajectory  can  stay  in  its  vicinity  for  only  a 
finite  amount  of  time  before  leaving  along  its  unstable  mani¬ 
fold  and  then,  enter  the  noisy  ball  at  the  original  periodic 
attractor  again,  and  so  on.  For  D^DC ,  the  probability  for  the 
trajectory  to  leave  the  noisy  ball  of  the  original  periodic  at¬ 
tractor  is  small.  Thus,  an  intermittent  behavior  can  be  ex¬ 
pected  where  the  trajectory  spends  long  stretches  of  time 
near  the  periodic  attractor,  with  occasional  bursts  out  of  it 
wandering  near  the  chaotic  saddle. 

A  consequence  of  the  noise-induced  intermittent  behav¬ 
ior  is  that  there  is  generally  unstable  dimension  variability 
associated  with  a  continuous  trajectory.  Under  noise,  both  the 
chaotic  saddle  and  the  original  periodic  attractor  belong  to  a 
single,  connected  dynamical  invariant  set.  Since,  in  the  ab¬ 
sence  of  noise,  periodic  orbits  on  the  chaotic  saddle  arc  all 
unstable  and  the  attractor  is  a  stable  periodic  orbit,  noise- 
induced  intermittency  means  that  a  trajectory  moves  in  re¬ 
gions  containing  periodic  orbits  with  distinct  unstable  dimen¬ 
sions.  A  feature  that  distinguishes  this  type  of  unstable 
dimension  variability  with  that  in  the  literature2-1  is  that 
here,  the  subsets  with  different  unstable  dimensions  are  lo¬ 
cated  in  distinct  regions  of  the  phase  space,  whereas  in  high¬ 
dimensional  chaotic  systems  such  as  the  kicked  double 
rotor,23'24  unstable  periodic  orbits  in  these  subsets  tend  to 


mix  with  each  other  densely  in  the  phase  space. 

At  a  fundamental  level,  the  appearance  of  unstable  di¬ 
mension  variability  implies  the  disappearance  of  the  neutral 
direction  of  the  flow.  Consider  a  three-dimensional  flow  in  a 
periodic  window,  where  the  periodic  attractor  contains  no 
unstable  direction  and  the  chaotic  saddle  possesses  one  un¬ 
stable  dimension.  The  role  of  noise,  when  it  is  sufficiently 
large  (£>>£>  t.),  is  to  link  these  two  dynamical  invariant  sets 
with  distinct  unstable  dimensions.  Now  examine  the  local 
eigenplanes  that  contain  the  neutral  direction  of  the  flow  as¬ 
sociated  with  the  periodic  attractor  and  the  chaotic  saddle.  In 
the  local  eigenplane  at  the  periodic  attractor,  there  is  a  stable 
direction  and  a  neutral  direction.  Consider  an  eigenvector  in 
the  neutral  direction.  In  the  eigenplane  of  a  point  in  the  cha¬ 
otic  saddle,  there  is  an  unstable  direction  and  a  neutral  direc¬ 
tion.  When  a  trajectory  is  driven  by  noise  from  the  periodic 
attractor  to  the  chaotic  saddle  along  its  stable  manifold,  the 
eigenvector  can  lie  anywhere  in  the  local  eigenplane  of  the 
corresponding  point  in  the  chaotic  saddle.  After  a  time,  the 
vector  will  be  aligned  in  the  unstable  direction,  due  to  the 
expanding  dynamics  of  the  chaotic  saddle.  Distances  along 
the  neutral  direction  of  the  original  periodic  attractor  can  no 
longer  be  preserved.  This  feature  of  a  noisy  chaotic  attractor 
is  fundamentally  different  from  that  of  a  deterministic  cha¬ 
otic  attractor,  where  a  neutral  direction  always  exists.  Thus 
we  see  that  unstable  dimension  variability  plays  a  fundamen¬ 
tal  role  in  shaping  the  topology  of  the  noisy  chaotic  flow. 

2.  Scaling  law  of  Lyapunov  exponents 

We  shall  argue  that  for  noise-induced  chaos,  the  largest 
Lyapunov  exponent  of  the  attractor  obeys  a  universal  alge¬ 
braic  scaling  law: 

XKDMD-D,)",  for  D^DC,  (14) 

where  the  scaling  exponent  a  depends  on  system  details. 
Consider  an  (N+  I  )-dimensional  flow  in  a  periodic  window. 
In  the  absence  of  noise,  the  chaotic  saddle  has  Ku  positive, 
one  zero,  and  Ks  negative  Lyapunov  exponents  (Ku  +  Ks 
=  N)  which  can  be  ordered  as  follows: 

=  \‘?()>  -  \f ■  > . . .  5*  -  \SK~- 1  s*  -  X'L  •  ( 1  5) 

The  periodic  attractor  has  one  zero  and  N  negative  expo¬ 
nents,  as  follows: 

0  =  \f’°>  —  \f~> - X,v“.  (16) 

For  D<DC,  an  asymptotic  trajectory  is  confined  in  the 
neighborhood  of  the  periodic  attractor,  so  the  largest 
Lyapunov  exponent  of  the  noisy  attractor  is  simply 
=  \/>0  =  0.  For  DsD,  (after  the  transition  to  chaos),  is 
approximately  given  by 

X^MD)\po+fs{D)\sKl=MD)\s/u  ,  (17) 

where  fP(D)  an d  fs(D)  are  the  probabilities  that  a  trajectory 
stays  near  the  original  periodic  attractor  and  the  chaotic 
saddle,  respectively.  Because  of  the  averaging  effect  of  noise, 
we  expect  the  dependence  on  noise  of  the  largest  Lyapunov 
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exponent  \£r  of  the  original  chaotic  saddle  to  be  weak.  Thus 
the  main  dependence  of  \\  on  noise  comes  from  /$(£>),  the 
frequency  of  visit  to  the  chaotic  saddle,  which  is  determined 
by  the  measure  of  its  stable  manifold  in  the  noisy  ball  of  the 
periodic  attractor. 

Consider  an  /V-dimensional  Poincare  map  corresponding 
to  the  (AT  1  )-dimensional  flow.  For  a  ball  of  radius  e ,  the 
natural  measure  of  the  stable  manifold  contained  within  is 
proportional  to 

^=(6Ay-v, 


where  eA'  is  proportional  to  the  volume  of  the  ball  and  ds  is 
the  information  dimension  of  the  stable  manifold  of  the  cha¬ 
otic  saddle.  Using  the  Kaplan  Yorke  conjecture21  for  chaotic 
saddles,4 1  ds  is  given  by 


ds=Ks+J+ 


//s-(\p+...+\5+) 

X  s  + 

KJ- M 


(18) 


where  J  is  an  integer  determined  by 


\r+...+xr+\y,= 


x 


and  Hs  is  the  forward  entropy  of  the  chaotic  saddle: 
1 

i-l  T 


(19) 


(20) 


Here  r  is  the  average  lifetime  of  the  chaotic  saddle  on  the 
Poincare  map.  (As  a  practical  matter,  r  is  in  the  unit  of  T,  the 
average  time  that  a  typical  trajectory  crosses  the  Poincare 
section.)  For  D<:DC ,  the  volume  of  the  noisy  ball  in  which 
the  stable  manifold  of  the  chaotic  saddle  lies  is  proportional 
to:  (Dn— £>a).  We  thus  have 

X  j  • —  (  Da  -  D;v)‘,«  /iV~  ( D  -  Dc)n . 


which  is  the  scaling  law  (14)  with  the  algebraic  scaling  ex¬ 
ponent  given  by 


Ks+J+ 


IIs—  (  X  ,  T  ...  T  X/^  ) 


(21) 


Numerical  support  for  the  scaling  law  can  be  found  in  Refs. 
31  and  42. 


B.  Noise  induced  bifurcation  and  chaos  in  mechanics 

We  continue  with  two  examples  of  how  small,  additive 
noise  induces  chaos  and  unstable  dimension  variability  in 
this  mechanical  system.  In  particular,  we  study  the  transition 
to  stochastic  chaos  when  we  add  stochastic  perturbations  of 
the  form 

^  =  F(y  ,p)  +  D£(t). 

where  F  is  the  deterministic  vector  field  of  the  mechanical 
system  defined  in  Eq.  (9),  p  represents  the  vector  of  param¬ 
eters,  and  D£(t)  is  the  additive  Gaussian  white  noise  with 
standard  deviation  D.  Note  that  £(/)  is  an  six-dimensional 
vector  whose  components  are  independent  Gaussian  random 
variables  of  zero  mean  and  scaled  variance.  Explicitly,  the 
noise  is  scaled  in  each  component  so  that  it  does  not  domi¬ 


nate  components  with  smaller  magnitudes.  We  find  the  ap¬ 
proximate  range  of  each  variable  when  it  is  in  its  steady  state 
and  then  scale  the  standard  deviation  of  the  noise  by  those 
factors  so  that  it  effects  each  component  equally.  Noise  is  not 
added  to  the  components  representing  the  drive  [^3  and 
in  Eq.  (9)],  so  those  variances  are  set  to  zero.  We  use  the 
standard  second-order  Milshtein  method,43  to  integrate  the 
stochastic  differential  equations.  We  also  integrate  the  Jaco¬ 
bian  and  calculate  the  Floquet  multipliers  to  find  the 
Lyapunov  exponents  based  on  time  averaging.  By  fixing  the 
parameters  of  Eq.  (9).  we  examine  the  dynamics  of  the  me¬ 
chanical  random  dynamical  system  as  the  noise  amplitude  is 
increased  from  zero. 

First,  we  study  the  system  with  small  coupling  using  the 
parameters  /x  =  0.086875  and  or  =1.8.  With  no  noise,  the 
dominant  behavior  of  the  system  is  a  chaotic  attractor  on  the 
slow  manifold  characterized  by  one  positive  Lyapunov  expo¬ 
nent.  As  expected  with  an  autonomous  flow,  the  second  larg¬ 
est  exponent  is  the  null  Lyapunov  exponent,  which  repre¬ 
sents  the  neutral  direction  associated  with  the  flow.  The  four 
remaining  Lyapunov  exponents  are  negative.  Adding  noise  to 
the  system  excites  a  nearby  high  dimensional  chaotic  saddle 
off  the  slow  manifold  and  the  potential  number  of  unstable 
(dynamical)  dimensions  is  increased  from  one  to  two.  Nu¬ 
merically,  we  observe  that  increasing  the  standard  deviation 
of  the  noise  (D)  increases  the  third  Lyapunov  exponent,  the 
largest  negative  exponent,  in  a  continuous  manner.  When  the 
third  exponent  is  close  to  crossing  zero,  the  null  exponent 
increases  away  from  zero,  leaving  no  zero  valued  Lyapunov 
exponent  within  a  small  window.  Define  the  beginning  of 
this  transition  Dc .  This  fundamentally  disturbs  the  noisy 
flow,  resulting  in  two  positive  Lyapunov  exponents  and  the 
third  largest  exponent  approaching  the  zero  value  from  the 
negative  side.  See  Fig.  13(a)  for  a  graph  of  this  transition. 
We  approximate  the  algebraic  scaling  of  the  Lyapunov  expo¬ 
nent  with  a  least  squares  fit  as  a~  1.5143,  having  a  maxi¬ 
mum  error  of  0.5786,  as  shown  in  Fig.  13(b). 

Due  to  the  small  coupling  parameter,  the  deterministic 
system  motion  is  constrained  to  a  slow  manifold,  which  is  a 
four-dimensional  surface.  The  chaotic  attractor  resides  on 
this  surface,  but  random  trajectories  experience  on  -off  inter- 
mittency,  which  is  characterized  by  a  bursting  behavior  off 
the  surface.  This  is  common  for  chaotic  attractors  having 
periodic  orbits  with  unstable  eigenvectors  transverse  to  the 
attractor.  As  an  orbit  approaches  the  attractor,  it  sometimes 
lands  near  one  of  these  repel lors,  which  ejects  it  from  the 
neighborhood  of  the  attractor.  Then,  the  trajectory  visits  the 
chaotic  saddle  off  the  slow  manifold  for  a  period  of  time 
until  it  begins  its  approach  back  to  the  attractor  once  again. 
Measuring  the  distance  of  a  trajectory  off  the  slow  manifold 
reveals  the  frequency  of  this  process,  and  an  average  bursting 
rate  can  be  calculated.  When  noise  is  added  to  the  system, 
the  bursting  rate  increases  with  the  standard  deviation  of  the 
noise.  This  is  expected  since  the  noise  facilitates  the  process 
of  a  trajectory  landing  near  a  repellor.  If  we  track  the  Euclid¬ 
ean  distance  of  the  trajectory  in  the  first  two  components 
from  the  slow  manifold,  we  can  set  a  threshold  to  define 
bursting  rates.  By  recording  the  fraction  of  iterates  in  a  long 
trajectory  as  a  function  of  the  standard  deviation  of  the  noise, 
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D  (standard  deviation) 


FIG.  13.  For  the  noisy  system  in  the  small  coupling  case,  (a)  the  second  and 
third  Lyapunov  exponents  versus  D  about  the  transition,  and  (b)  algebraic 
scaling  of  the  second  largest  Lyapunov  exponent  with  D  —  Dc.  The  solid 
line  in  indicates  the  theoretical  slope  ci. 

we  see  a  change  in  behavior  at  the  same  critical  value  Dc , 
where  the  number  of  positive  Lyapunov  exponents  increase. 
See  Fig.  14.  The  data  follow  exponential  growth,  and  using  a 
least  squares  fit.  we  can  fit  line  to  the  natural  log  of  the  data 
with  slope  1 1.6296,  which  has  a  maximum  error  of  0.2148, 
as  shown  in  Fig.  14. 

Another  quantity  that  we  can  measure  is  the  interspike 
interval.  This  is  the  time  between  bursts  away  from  the  slow 
manifold.  Recording  the  interspike  intervals  along  a  random 
trajectory  also  indicates  how  noise  increases  the  frequency  of 
transients  on  the  chaotic  saddle.  We  calculate  the  interspike 
interval  using  the  same  threshold  as  the  burst  rate.  Because 
the  distance  measurement  does  not  have  a  pattern  in  time  or 
amplitude,  it  is  difficult  to  define  the  beginning  and  end  of  an 
individual  burst.  Therefore,  we  define  the  beginning  of  each 
burst  as  the  time  when  the  distance  increases  above  the  burst 
threshold.  The  interspike  intervals  are  the  intervals  between 
these  times.  Consider  the  sequence  of  interspike  intervals  for 
a  random  trajectory  using  the  parameters  fi~  0.086  875,  a 
=  1.8,  and  D^0.25.  The  histogram  of  interspike  interval 
(1SI)  sequence  follows  a  semilog  scaling  law.  The  slope  of 
the  least  squares  fit  is  -0.1275  with  a  maximum  error  of 
1.9619,  as  shown  in  Fig.  15. 

At  larger  coupling  parameters,  the  dominant  behavior  of 
the  system  is  a  periodic  attractor.  Therefore,  with  no  noise, 
all  the  Lyapunov  exponents  are  negative  except  for  the  null 
exponent.  The  addition  of  noise  emulates  chaotic  behavior, 
which  is  observed  by  the  bifurcation  in  the  Lyapunov  expo¬ 
nents  from  zero  to  two  local  unstable  directions.  We  use  the 
parameters  /x  =  0.5025  and  a-^0.65  as  an  example  in  Fig. 
16.  Notice  that  the  two  largest  exponents  increase  in  a  con¬ 


FIG.  14.  For  the  noisy  system  using  the  parameters  /a  =  0.086  875  and  a 
=  1.8,  (a)  the  bursting  rate,  and  (b)  algebraic  scaling  of  the  bursting  rate. 
The  solid  line  in  indicates  the  slope. 


tinuous  manner  above  zero  and  the  third  largest  approaches 
zero  from  the  negative  side  in  Fig.  17(a).  This  is  a  new  type 
of  transition  to  noise  induced  chaos.  The  largest  Lyapunov 
exponent  follows  an  algebraic  scaling  similar  to  the  one 
shown  in  the  small  coupling  case.  We  use  Dc—  0.002  938, 
which  results  in  a  linear  fit  with  slope  1,4973  and  maximum 
error  of  0.8087.  The  second  largest  exponent  also  follows  an 
algebraic  scaling,  but  it  increases  from  a  negative  value. 
Therefore,  we  must  translate  the  Lyapunov  exponent  by  that 
negative  value  so  we  can  use  the  natural  log.  We  use  the 
value  of  the  Lyapunov  exponent  at  Dc ,  called  Lc ,  and  find 
ln(L-Lc.)  as  we  increase  D  from  Dc .  This  results  in  a  linear 
fit  with  slope  1.5090  and  maximum  error  of  0.8198.  Note  the 
similarity'  in  this  scaling  to  that  of  the  maximum  Lyapunov 
exponent.  Graphs  of  each  of  these  fits  are  shown  in  Fig. 


IS1  histogram 


FIG.  15.  The  interspike  interval  statistics  for  the  parameters  fx 
=  0.086  875,  «=1.8,  and  0  =  0.25. 
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FIG.  16.  The  Lyapunov  exponents  of  the  noisy  system.  This  is  the  large 
coupling  case  using  the  parameters  fx~ 0.5025  and  or  — 0.65. 


17(b).  Because  this  example  uses  a  large  coupling  parameter, 
the  dynamics  are  not  dominated  by  the  slow  manifold  and 
we  cannot  compute  bursting  statistics  similar  to  the  previous 
example. 

V.  CONCLUSIONS 

We  have  explored  a  class  of  linear  continuum  mechani¬ 
cal  systems  coupled  to  a  nonlinear  oscillator  from  a  dynami¬ 
cal  systems  point  of  view.  In  previous  work,  such  systems 
were  ideally  driven  from  an  outside  source.  In  the  model 
presented  here,  we  considered  a  more  general  case  where  the 
frequency  of  the  driver  is  slightly  perturbed  by  the  momen¬ 
tum  change  of  the  mechanical  system  to  which  it  is  con¬ 
nected.  We  have  derived  an  ODE-PDE  system  describing 
the  physics,  and  then  showed  how  to  construct  the  modal 
decomposition  into  an  infinite  set  of  ODE’s.  A  truncated 


ln(D-Dc) 


FIG.  17.  These  graphs  use  the  noisy  system  with  large  coupling  parameters: 

=  0.5025  and  a  =  0.65.  (a)  A  close  up  of  the  three  largest  Lyapunov  ex¬ 
ponents  near  the  transition  to  stochastic  chaos,  (b)  The  algebraic  scaling  for 
the  two  largest  Lyapunov  exponents.  The  largest  exponent  is  shown  in  black 
and  the  data  arc  plotted  ln(/_)  vs  ln(D-Df  ).  The  second  largest  is  in  gray  and 
the  data  arc  plotted  In (L~L(.)  vs  In Both  linear  fits  have  slope  close 
to  1.5. 


model  was  then  derived  for  analysis,  in  which  the  system 
was  described  as  a  slow  •  fast  time  scale  system,  which  could 
be  analyzed  using  singular  perturbation  theory. 

The  deterministic  system  was  analyzed  for  its  bifurca¬ 
tion  structure,  which  included  routes  to  chaos,  as  well  as 
torus  bifurcations.  However,  one  of  the  interesting  bifurca¬ 
tions  observed  was  that  of  a  dimension  changing  bifurcation. 
When  the  singular  parameter  (fi)  was  small,  one  could  ob¬ 
serve  chaos  constrained  in  a  neighbor  of  an  underlying  in¬ 
variant  manifold.  However,  increasing  a  parameter,  such  as 
amplitude  a  drive  term  in  Eq.  (9),  would  cause  the  dynamics 
to  burst  off  the  manifold  into  the  ambient  space.  Accompa¬ 
nying  sufficient  bursting  was  the  appearance  of  an  additional 
Lyapunov  exponent,  signaling  a  change  in  dynamical  dimen¬ 
sion. 

Since  the  invariant  slow  manifold  is  an  unstable  object  at 
appropriate  parameters,  constraining  the  dynamics  to  the 
manifold  reveals  a  chaotic  saddle  structure.  Such  a  structure 
is  important  when  examining  both  deterministic  and  stochas¬ 
tic  bifurcations  from  low  to  high  dimensions. 

One  important  aspect  of  the  noise  induced  bifurcation 
was  quantifying  the  relationship  of  the  change  in  Lyapunov 
exponents  with  respect  to  the  standard  deviation  of  the  noise. 
Here  we  used  the  unstable  dimension  variability  of  the  sys¬ 
tem  to  show  explicitly  in  a  mechanical  system  that  the  expo¬ 
nent  obeys  a  universal  scaling  law.  It  is  a  remarkable  fact  that 
the  scaling  actually  persists  over  a  wide  range  of  noise  am¬ 
plitudes.  It  is  also  an  interesting  observation  that  the  change 
is  continuous,  even  when  the  deterministic  changes  from  pe¬ 
riodic  to  hyper-chaotic  behavior  are  discontinuous. 

Several  areas  of  inquiry  are  suggested  by  the  current 
study.  First,  recent  work  has  shown  that  sudden  changes  may 
occur  as  a  dimension  changing  bifurcation  when  mechanical 
systems  of  sufficient  complexity  are  operated  near  resonance. 
Since  we  observe  such  changes  in  simple  prototypical  sys¬ 
tems,  as  we  do  here,  it  appears  that  many  other  systems  may 
exhibit  the  same  nonlinear  vibrations.  If  so,  novel  controls 
may  be  designed  based  on  the  invariant  manifold  theory. 
Moreover,  such  control  may  be  used  to  spread  energy  into 
higher  heat  dissipative  modes,  or  other  governing  devices, 
making  energy  transfer  from  one  part  of  a  structure  to  an¬ 
other  more  efficient  and  directed. 

On  the  other  hand,  since  additive  noise  modifies  the  dy¬ 
namics  and  its  dimension  continuously,  it  may  be  used  a 
better  probe  for  nondestructive  evaluation.  By  comparing  the 
pristine  system  to  later  use.  one  may  use  time  series  tests  for 
nonstationaritv  to  explore  controlled  frequency  responses  by 
adding  noise. 

Other  areas  of  interest  include  nonlinear  elasticity,  dif¬ 
ferent  boundary  conditions,  such  as  clamping,  and  higher 
dimensional  structures,  such  as  trusses. 
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